Possibility of direct Mott insulator-to-superfiuid transitions 
in weakly disordered boson systems 
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We study the zero-temperature phase transitions of a two-dimensional disordered boson Hubbard 
model at incommensurate boson densities. Via matrix diagonalization and quantum Monte Carlo 
simulations, we construct the phase diagram and evaluate the correlation length exponent v. In the 
presence of weak disorder, we obtain v — 0.5 ± 0.1, the same value as that in the pure model, near 
the tip of a Mott insulator lobe, using the dynamical critical exponent 2 = 2. As the strength of 
disorder is increased beyond a certain value, however, the value of v is found to change to 0.9 ± 0.1. 
This result strongly suggests that there exist direct Mott insulator-to-superfluid transitions around 
the tip of a Mott insulator lobe in the weak disorder regime. 
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Two-dimensional (2D) interacting boson systems dis- 
play quantum phase transitions from an insulating 
state to a superfluid (SF) state at zero temperature. 
Physical realization of this transition may include dis- 
ordered thin film superconductors, Josephson-junction 
arrays, granular superconductors and 4 He films ad- 
sorbed in porous media In the absence of disorder 
the insulating state is a Mott insulator (MI), which has 
a commensurate value of the boson density since, other- 
wise, excessive bosons or holes will move freely to yield 
a superfluid. The excessive bosons or holes can be lo- 
calized, on the other hand, in the presence of disorder, 
to make an incommensurate insulator. The resulting 
Bose glass (BG) insulating phase has attracted consid- 
erable interest recently. In the mean-field theory jjj, 
it has been argued that in the presence of disorder the 
transition from a Mott insulator to a superfluid occurs 
only through the BG phase. However, a recent quantum 
Monte Carlo study ||, performed at the tip of an MI 
lobe (i.e., at a commensurate value of the boson density), 
has demonstrated the existence of a direct MI-SF tran- 
sition. The subsequent renormalization group study || 
has been interpreted to suggest that the direct MI-SF 
transition occurs around the tip of the MI lobe in high di- 
mensions (d > 4) but only at the tip in lower dimensions 
(2 < d < 4). These results raise an interesting question 
as to whether such a direct MI-SF transition occurs even 
at incommensurate values of the boson density. 

This work investigates the possibility of the direct MI- 
SF transition off the tip of an MI lobe in two dimensions. 
We perform both matrix diagonalization and quantum 
Monte Carlo simulations, and find evidences for the di- 
rect transition around the tip of an MI lobe in the weak 
disorder limit. The superfluid onset points are identi- 
fied and the correlation length critical exponents are es- 
timated at various disorder strengths. The corresponding 
phase diagram for weak disorder is constructed. 

We consider the boson Hubbard model with disorder, 



described by the Hamiltonian 
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where b\ and bi are the boson creation and destruc- 
tion operators at site i on an L x L square lattice, and 
rii = b\bi is the number operator. In Eq. (|l|), U is 
the strength of the on-site repulsion, /i is the chemical 
potential, Vi is the random on-site potential distributed 
uniformly between —A and A, and t measures the hop- 
ping strength between nearest neighboring sites. In the 
limit of the large number of particles, we may take the 
phase-only approximation [Q and reduce Eq. ([!]) to the 
quantum phase Hamiltonian 
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where fa is the phase of the bosons condensed at site i 
and satisfies the relation [rij, fa] — idij. Note that in this 
representation, rij denotes the deviation from the mean 
integer number no (around the mean number (b\bi)). 

In order to determine the phase boundary, we first use 
the matrix diagonalization method with a truncated ba- 
sis set. The zero-temperature phase boundary between 
the MI phase and the SF/BG phase is determined by 
comparing the ground state energy of the system at com- 
mensurate boson density and that of the system contain- 
ing an extra hole or particle. The basis set is chosen 
to include the lowest-energy states of the Hamiltonian 
given by the first and the second term in Eq. (||) (i.e., 
the zeroth order in t) and the states coupled to these by 
hopping up to (2n+l)th order in t, where we set n = 1 
in this work. We identify the lowest-energy state among 
the states of which the total boson number is zero as the 
MI state and the lowest one among the states with the 
total boson number being unity as the SF/BG state, and 
determine the phase boundary by locating the values of 



fi (> 0) and t at which the energies of the two states 
become equal. 

Since the basis set includes those states coupled up to 
the 3rd order in t, the result will be the same as that from 
the energy perturbation to the same order. The numbers 
of states involved in this calculation are 4L 4 — 9L 2 and 
4L 2 + 1 in the cases of the BG/SF state and of the MI 
state, respectively. We adopt the Lanczos method [jy] to 
obtain the lowest eigenvalue and the corresponding eigen- 
state, and take 2000 disorder realizations for each system 
size. In the pure case (A = 0), we find the phase diagram 
which is consistent with the previous perturbation result 
of Freericks and Monien |Q . 

Figure [l] shows the zero temperature phase diagram 
with the disorder strength A = 0.2 up to the system 
size L = 10 obtained from the truncated basis set. The 
shape of the MI lobe constructed from the previous per- 
turbation result in the pure case [ fl0| should be round 
and shorter than that obtained here. In the absence of 
hopping (t = 0), where no SF state is possible, the phase 
boundary (between the MI and BG states) approaches, 
as the system size L is increased, 0.5— A, as it should. 

To distinguish the SF phase from the BG phase, we 
use localization argument, and define the participation 
ratio as 



PL 
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where pt is the probability that particles are found at site 
i ||, and [• • -] av denotes the average over different dis- 
order realizations. At the generic SF-BG transition, the 
participation ratio satisfies the scaling relation 



(4) 



where 5 t — (t — t c )/t c is the distance from the critical 
point t c and v is the correlation length exponent. Here 
the scaling function p and an additional exponent y have 
been introduced jfjjl^]. Note that since we just consider 
the states in which the total number of bosons is unity, 
this method is useful to determine the critical point very 
close to the MI phase, i.e., on the phase boundary in 

Fig '0- n 

Figure presents the finite-size scaling of the partici- 
pation ratio of the SF/BG state. From this scaling be- 
havior, we obtain y = 0.92 ± 0.03, v = 1.5 ± 0.3, and 
t c = 0.12 ± 0.02; the latter separates the SF state from 
the BG one as t (or fi) is varied along the phase boundary 
in Fig. |l[ Accordingly, the system undergoes a direct MI- 
SF transition as the phase boundary is crossed for small 
/i (< 0.18 in Fig. [l]), i.e., near the tip of an MI lobe. Note 
that the possible origin of this direct MI-SF transition, as 
discussed later, is the abundance of particle-hole excita- 
tions around the tip of a lobe. This suggests that the size 
of the basis set including the states up to the 3rd order in 
t might not be sufficient to determine the phase bound- 
ary unambiguously and to compute the corresponding 



exponents accurately. Such abundance of particle-hole 
excitations is expected to forbid the perturbation calcu- 
lation to be accurate, possibly causing the discrepancy 
between the perturbative results ||Io|] and the quantum 
Monte Carlo results || in the pure case without disor- 
der. 

It is thus needed to investigate the transition via quan- 
tum Monte Carlo simulations, which are in general more 
reliable. For that purpose, we follow the standard proce- 
dure ||[)| to transform the 2D quantum phase Hamilto- 
nian in Eq. ^) to the (2+l)-dimensional classical action 
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integer current vector J 



(i,t) 



where the 

(J(i t)' J(i t)' J(i t)) 1S divergenceless on each lattice site 
(i,t) as indicated. The coupling constant K, correspond- 
ing roughly to y/t/U, takes the role of the temperature. 

We perform quantum Monte Carlo simulations on the 
classical action in Eq. (|B|), employing the heat bath al- 
gorithm [|| at a classical temperature K. An important 
quantity in the analysis is the zero-frequency superfluid 
stiffness 
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where n, 



(V-^) 2(i t) J?i t) ^ s t ne winding number 
along the (spatial) x direction. The finite-size scaling 
behavior of the superfluid stiffness P,[L3| reads 
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where 5 = (K—K c )/K c is the distance from the critical 
point K Cl the spatial dimension d is two in this work, and 
z is the dynamical critical exponent. 

In order to investigate the scaling behavior, the value 
of z should be known in advance. At the very tip of an MI 
lobe, the direct MI-SF transition in the presence of weak 
disorder shows the same behavior as the pure system: the 
dynamical critical exponent z = 1 and the correlation ex- 
ponent v = 0.67 ||. Off the tip, we expect the dynamical 
exponent z = 2 in the disordered system since the com- 
pressibility is finite at the transition point [Q; the same 
number is expected even in the pure case. We thus set 
the dynamical exponent z — 2, and measure the correla- 
tion length exponent i/, the value of which was estimated 
in the previous studies to be 0.9±0.1 in the BG-SF tran- 
sition HH and 0.5 ± 0.1 in the MI-SF transition 

Keeping L T /L Z constant, we simulate the systems of 
sizes LxLxL T = 6x6x9, 8x8x16, and 10 x 10 x 25. To 
tune the transition, we vary the temperature K while fix- 
ing /J, and A. We further take the average over 60 — 200 
disorder realizations and, for each disorder realization, 
perform typically 4000 — 80000 Monte Carlo sweeps for 
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equilibration, followed by equally many sweeps for mea- 
surement. The equilibration is checked through the use 
of the standard equilibration test technique |9||l4|] . 

In Fig. ||, we show the results with z = 2 at fi = 0.3 and 
A = 0.1. One can identify clearly the common crossing 
point of L z p at K c = 0.257 ± 0.003 as the critical point. 
The inset of Fig. || shows a scaling plot, which yields 
v = 0.5 ±0.1, the same value as the pure model Q. On 
the other hand, in the strong-disorder case (A = 0.3) 
shown in Fig. |[ the scaling behavior yields v — 0.9±0.1, 
which agrees well with the previous results of the BG-SF 
transition |^|J|,[|. These demonstrate that in the weak 
disorder regime a direct MI-SF transition takes place not 
only at the tip but also off the tip of the MI lobe. Figure 
|H| presents the critical temperature K c as a function of 
the disorder strength A for fj, = 0.3, displaying the direct 
MI-SF transition occurring for A < 0.16. It is of inter- 
est that the point on which the value of the correlation 
length exponent v changes appears to have the maximum 
slope on the K c -A curve. 

Figure ^ summarizes the phase diagram of the dis- 
ordered boson Hubbard model at the disorder strength 
A = 0.2. The phase boundary between the MI and the 
BG phases is estimated by means of the matrix diago- 
nalization in the system of size L — 10. Here we have 
used the results in Ref. J|,[l(| to scale K as a function of 
t, which is valid only for small t. 

Finally, we discuss the possible origin of the direct 
MI-SF transition around the tip of an MI lobe. One 
might think that off the tip the BG phase intervenes 
very slimly between the MI and SF phases and that our 
results supporting the direct MI-SF transition merely re- 
flect finite-size effects. Be this the case, the localization 
length would be large near the tip and the BG phase 
would manifest itself only on the length scale exceeding 
the system size. Then, as pointed out in Ref. ||, an 
anomalous behavior is expected to occur to change the 
correlation exponent v, suggesting continuous change of 
the exponent with the thickness of the BG phase near 
the tip. On the other hand, our results in Figs. |5| and ^ 
indicate that the exponent changes quite abruptly, mak- 
ing the above scenario rather unlikely. Near the tip, the 
particle-hole excitations are ubiquitous and presumably 
tend to suppress the disorder effects, allowing the pos- 
sibility of an extended state for the extra boson. Thus 
the resulting direct MI-SF transition around the tip may 
reflect the peculiar nature of boson localization. 

In summary, we have studied the disordered boson 
Hubbard model by means of the matrix diagonalization 
with restricted basis states, which include those states 
overlapping with each other through nearest-neighbor 
hopping up to the 3rd order. The finite-size scaling of 
the participation ratio at A = 0.2 gives evidence for the 
direct MI-SF transition at an incommensurate value of 
the boson density. To investigate this direct transition 
more clearly, we have investigated the scaling behavior of 
the superfluid stiffness via quantum Monte Carlo simu- 
lations and found that as the disorder strength is varied, 



the value of the correlation exponent v changes rather 
abruptly from the weak-disorder value 0.5 to the strong- 
disorder one 0.9. This indicates that the direct MI-SF 
transition occurs for weak disorder at an incommensu- 
rate density. The possible origin of this transition could 
be the abundance of the particle-hole excitations around 
the tip of an MI lobe, which suggests the peculiar nature 
of boson localization. 
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FIG. 1. Phase diagram in the (t, n) plane for the disorder 
strength A = 0.2. Each line gives an estimate of the phase 
boundary between the MI state and the SF or BG state, as 
determined by matrix diagonalization. 
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FIG. 4. Behavior of the stiffness p for fj, = 0.3 and 
A = 0.3, with the aspect ratio kept constant for z = 2. The 
inset shows a scaling plot with K c = 0.223 and v — 0.9. 




FIG. 2. Scaling plot of the participation ratio p L near the FIG. 5. Critical point K c as a function of the disorder 

MI-SF/BG phase boundary, which yields t c = 0.12, y = 0.92, strength A for fi = 0.3. The solid dots indicate the direct 
and v = 1.5. MI-SF transition points and the open dots the BG-SF tran- 

sition points. 




A = 0.1, with the aspect ratio kept constant for 2 = 2. The K 
inset shows a scaling plot with K c = 0.257 and v = 0.5. 
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FIG. 6. Phase diagram in the (K,/j,) plane for A = 0.2. 
The dotted line (without symbols) corresponds to the rescaled 
data from the matrix diagonalization method, while the solid 
and dashed lines (with symbols) are the results of quantum 
Monte Carlo simulations. The two points marked by crosses 
represent the results in Ref. [5]. 



